home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / graphic / pvquan16.zip / QUANT / HECKBERT.C < prev    next >
C/C++ Source or Header  |  1992-11-30  |  17KB  |  568 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  *                  Copyright (c) 1991, Frank van der Hulst             *
  4.  *                          All Rights Reserved                         *
  5.  *                                                                      *
  6.  * Authors:                                                             *
  7.  *        FvdH - Frank van der Hulst (Wellington, NZ)                   *
  8.  *                                                                      *
  9.  * Versions:                                                            *
  10.  *    V1.1 910626 FvdH - QUANT released for DBW_RENDER                  *
  11.  *    V1.2 911021 FvdH - QUANT released for PoV Ray                     *
  12.  *    V1.4 920303 FvdH - Ported to GNU                                  *
  13.  *    V1.6 921023 FvdH - Produce multi-image GIFs                       *
  14.  *                     - Port to OS/2 IBM C Set/2                       *
  15.  *                                                                      *
  16.  ************************************************************************/
  17. /*
  18.  * This software is copyrighted as noted below.  It may be freely copied,
  19.  * modified, and redistributed, provided that the copyright notice is
  20.  * preserved on all copies.
  21.  *
  22.  * There is no warranty or other guarantee of fitness for this software,
  23.  * it is provided solely "as is".  Bug reports or fixes may be sent
  24.  * to the author, who may or may not act on them as he desires.
  25.  *
  26.  * You may not include this software in a program or other software product
  27.  * without supplying the source, or without informing the end-user that the
  28.  * source is available for no extra charge.
  29.  *
  30.  * If you modify this software, you should include a notice giving the
  31.  * name of the person performing the modification, the date of modification,
  32.  * and the reason for such modification.
  33. */
  34. /*
  35.  * colorquant.c
  36.  *
  37.  * Perform variance-based color quantization on a "full color" image.
  38.  * Author:    Craig Kolb
  39.  *        Department of Mathematics
  40.  *        Yale University
  41.  *        kolb@yale.edu
  42.  * Date:    Tue Aug 22 1989
  43.  * Copyright (C) 1989 Craig E. Kolb
  44.  * $Id: colorquant.c,v 1.3 89/12/03 18:27:16 craig Exp $
  45.  *
  46.  * $Log:    colorquant.c,v $
  47.  *
  48.  * Revision 1.4  91/06/26  16:00:00  Frank van der Hulst
  49.  * Ported to Turbo C;
  50.  *       Call farmalloc rather than malloc
  51.  *       Virtual memory added to swap box data to/from disk
  52.  *       Rewritten in ANSI C
  53.  *       Removed call to QuantHistogram() from colorquant, to allow two
  54.  *       image files to create one palette
  55.  *       Changed QuantHistogram() to read from file, rather than from an
  56.  *       array
  57.  *          Changed format of palette to conform with the VGA palette
  58.  *
  59.  * Revision 1.3  89/12/03  18:27:16  craig
  60.  * Removed bogus integer casts in distance calculation in makenearest().
  61.  *
  62.  * Revision 1.2  89/12/03  18:13:12  craig
  63.  * FindCutpoint now returns FALSE if the given box cannot be cut.  This
  64.  * to avoid overflow problems in CutBox.
  65.  * "whichbox" in GreatestVariance() is now initialized to 0.
  66.  *
  67.  */
  68.  
  69. #ifdef __TURBOC__
  70. #include <mem.h>
  71. #define HUGE 1.79e308
  72. #endif
  73. #include <math.h>
  74.  
  75. #include "quant.h"
  76. #include "heckbert.h"
  77.  
  78. #define MAX(x,y)    ((x) > (y) ? (x) : (y))
  79.  
  80. static unsigned long        NPixels = 0L;            /* total # of pixels */
  81.  
  82. static int neighbours[MAXCOLORS];
  83.  
  84. /*
  85.  * Compute the histogram of the image as well as the projected frequency
  86.  * arrays for the first world-encompassing box.
  87.  * We compute both the histogram and the proj. frequencies of
  88.  * the first box at the same time to save a pass through the
  89.  * entire image.
  90.  * The projected frequency arrays of the largest box are zeroed out as
  91.  * as part of open_box_file(), called from main().
  92.  */
  93.  
  94. void QuantHistogram(Box *box)
  95. {
  96. unsigned long *rf, *gf, *bf;
  97. UCHAR pixel[3];
  98.  
  99.     rf = box->freq[RED];
  100.     gf = box->freq[GREEN];
  101.     bf = box->freq[BLUE];
  102.  
  103.     while (get_pixel(pixel)) {
  104.         rf[pixel[RED]]++;
  105.         gf[pixel[GREEN]]++;
  106.         bf[pixel[BLUE]]++;
  107.         Histogram[(((pixel[RED]<<INPUT_BITS)|pixel[GREEN])<<INPUT_BITS)|pixel[BLUE]]++;
  108.         NPixels++;
  109.     }
  110. }
  111.  
  112. /*
  113.  * Compute mean and weighted variance of the given box.
  114.  */
  115. void BoxStats(Box HUGE_PTR box)
  116. {
  117. int i, color;
  118. unsigned long *freq;
  119. double mean, var;
  120.  
  121.     if(box->weight == 0) {
  122.         box->weightedvar = 0;
  123.         return;
  124.     }
  125.  
  126.     box->weightedvar = 0.0;
  127.     for (color = 0; color < 3; color++) {
  128.         var = mean = 0;
  129.         i = box->low[color];
  130.         freq = &box->freq[color][i];
  131.         for (; i < box->high[color]; i++, freq++) {
  132.             mean += i * *freq;
  133.             var += i*i* *freq;
  134.         }
  135.         box->mean[color] = mean / box->weight;
  136.         box->weightedvar += var - box->mean[color]*box->mean[color]*box->weight;
  137.     }
  138.     box->weightedvar /= NPixels;
  139. }
  140.  
  141. /*
  142.  * Return the number of the box in 'boxes' with the greatest variance.
  143.  * Restrict the search to those boxes with indices between 0 and n-1.
  144.  */
  145. int GreatestVariance(int n)
  146. {
  147.     int i, whichbox = 0;
  148.     double max;
  149.     Box *box;
  150.  
  151.     max = -1;
  152.     for (i = 0; i < n; i++) {
  153.         box = get_box_tmp(i);
  154.         if (box->weightedvar > max) {
  155.             max = box->weightedvar;
  156.             whichbox = i;
  157.         }
  158.     }
  159.     return whichbox;
  160. }
  161.  
  162. /*
  163.  * Update projected frequency arrays for two boxes which used to be
  164.  * a single box.
  165.  */
  166. void UpdateFrequencies(Box HUGE_PTR box1, Box HUGE_PTR box2)
  167. {
  168. unsigned long myfreq, HUGE_PTR h;
  169. int b, g, r;
  170. int roff;
  171.  
  172.     memset(box1->freq[0], 0, IN_COLOURS * sizeof(unsigned long));
  173.     memset(box1->freq[1], 0, IN_COLOURS * sizeof(unsigned long));
  174.     memset(box1->freq[2], 0, IN_COLOURS * sizeof(unsigned long));
  175.  
  176.     for (r = box1->low[0]; r < box1->high[0]; r++) {
  177.         roff = r << INPUT_BITS;
  178.         for (g = box1->low[1];g < box1->high[1]; g++) {
  179.             b = box1->low[2];
  180.             h = Histogram + (((roff | g) << INPUT_BITS) | b);
  181.             for (; b < box1->high[2]; b++) {
  182.                 if ((myfreq = *h++) == 0)
  183.                     continue;
  184.                 box1->freq[0][r] += myfreq;
  185.                 box1->freq[1][g] += myfreq;
  186.                 box1->freq[2][b] += myfreq;
  187.                 box2->freq[0][r] -= myfreq;
  188.                 box2->freq[1][g] -= myfreq;
  189.                 box2->freq[2][b] -= myfreq;
  190.             }
  191.         }
  192.     }
  193. }
  194.  
  195. /*
  196.  * Compute the 'optimal' cutpoint for the given box along the axis
  197.  * indicated by 'color'.  Store the boxes which result from the cut
  198.  * in newbox1 and newbox2.
  199.  */
  200. int FindCutpoint(Box HUGE_PTR box, int color, Box HUGE_PTR newbox1, Box HUGE_PTR newbox2)
  201. {
  202.     double u, v, max;
  203.     int i, maxindex, minindex, cutpoint;
  204.     unsigned long optweight, curweight;
  205.  
  206.     if (box->low[color] + 1 == box->high[color])
  207.         return FALSE;    /* Cannot be cut. */
  208.     minindex = (int)((box->low[color] + box->mean[color]) * 0.5);
  209.     maxindex = (int)((box->mean[color] + box->high[color]) * 0.5);
  210.  
  211.     cutpoint = minindex;
  212.     optweight = box->weight;
  213.  
  214.     curweight = 0.;
  215.     for (i = box->low[color] ; i < minindex ; i++)
  216.         curweight += box->freq[color][i];
  217.     u = 0.;
  218.     max = -1;
  219.     for (i = minindex; i <= maxindex ; i++) {
  220.         curweight += box->freq[color][i];
  221.         if (curweight == box->weight)
  222.             break;
  223.         u += (double)(i * box->freq[color][i]) / box->weight;
  224.         v = ((double)curweight / (box->weight-curweight)) *
  225.                 (box->mean[color]-u)*(box->mean[color]-u);
  226.         if (v > max) {
  227.             max = v;
  228.             cutpoint = i;
  229.             optweight = curweight;
  230.         }
  231.     }
  232.     cutpoint++;
  233.     *newbox1 = *newbox2 = *box;
  234.     newbox1->weight = optweight;
  235.     newbox2->weight -= optweight;
  236.     newbox1->high[color] = cutpoint;
  237.     newbox2->low[color] = cutpoint;
  238.     UpdateFrequencies(newbox1, newbox2);
  239.     BoxStats(newbox1);
  240.     BoxStats(newbox2);
  241.  
  242.     return TRUE;    /* Found cutpoint. */
  243. }
  244.  
  245. /*
  246.  * Cut the given box.  Returns TRUE if the box could be cut, FALSE otherwise.
  247.  */
  248. int CutBox(Box HUGE_PTR box, Box HUGE_PTR newbox)
  249. {
  250.     int i;
  251.     double totalvar[3];
  252.     static Box newboxes[3][2];  /* Only used by CutBox, but don't want it on stack */
  253.  
  254.     if (box->weightedvar == 0. || box->weight == 0)
  255.         /*
  256.          * Can't cut this box.
  257.          */
  258.         return FALSE;
  259.  
  260.     /*
  261.      * Find 'optimal' cutpoint along each of the red, green and blue
  262.      * axes.  Sum the variances of the two boxes which would result
  263.      * by making each cut and store the resultant boxes for
  264.      * (pos